% Solve for second best allocation
%
% Step 0   : Guess a Lagrange multiplier xi for the RC
% Step 1   : Solve for the optimal allocation
%      1-1 : Guess c1
%      1-2 : Compute the allocation 
%      1-3 : Check the IC for the highest type. Iterate until convergence
% Step 2   : Compute the residual of the RC. Iterate until convergence 

temp.xi  = fzero(@(x) SecondBestStep2(x,pr,pw,w,pai,G),1.4038708,optimset('TolX',1e-9));
SecondBestStep2(temp.xi,pr,pw,w,pai,G); 

%% Allocation
load c
load y
sb.c       = c;
sb.y       = y;
sb.output  = sum(pai.*sb.y);
sb.welfare = sum(pw.*pai.*Utility(sb.c,sb.y,w,pr));
sb.MTR     = 1-pr.Om./sb.c.^(-pr.gam).*w.^-(1+pr.sig).*sb.y.^pr.sig;
sb.ATR     = (sb.y-sb.c)./sb.y;
sb.D       = zeros(pr.na,1);
for i = 1:pr.na
    sb.D(i)= 1-sum(pai(i:pr.na).*pw(i:pr.na).*sb.c(i:pr.na).^-pr.gam)./sum(pai(i:pr.na))./sum(pw.*sb.c.^-pr.gam.*pai);
end
sb.Dhazard = zeros(pr.na,1);
for i = 1:pr.na
    sb.Dhazard(i)= sb.D(i).*sum(pai(i:pr.na))./pdf(i);
end
sb.WelfareGain = 100*fzero(@(x) sum(pw.*pai.*Utility((1+x).*baseline.c,baseline.y,w,pr))-sb.welfare,.5);
sb.WelfareGain = 100*fzero(@(x) sum(pw.*pai.*Utility((1+x).*baseline.c,baseline.y,w,pr))-sb.welfare,.5);

sb.OutputLoss  = 100*(sb.output/baseline.output-1);
sb.EfAveMTR    = sum(pai.*sb.MTR.*sb.y/sb.output);
[temp.y,temp.i]= min(sb.y);
temp.c         = sb.c(temp.i);
sb.Tr2Y        = (temp.c-temp.y)./sb.output;
sb.TrG2Y       = (temp.c-temp.y+G)./sb.output;

save ./variables/sb sb
clear temp; delete c.mat y.mat